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Abstract.  We  propose  to  perform  turbulent  flow  simulations  in  such  a 
manner  that  the  difference  operators  do  have  the  same  symmetry  prop- 
erties as  the  underlying  differential  operators,  i.e.  the  convective  operator 
is  represented  by  a skew-symmetric  matrix  and  the  diffusive  operator  is 
approximated  by  a symmetric,  positive-definite  matrix.  Such  a symmetry- 
preserving discretization  of  the  Navier-Stokes  equations  is  stable  on  any 
grid,  and  conserves  the  total  mass,  momentum  and  kinetic  energy  (when 
the  physical  dissipation  is  turned  off).  Its  accuracy  is  tested  for  a turbu- 
lent channel  flow  at  Re=5,600  (based  on  the  channel  width  and  the  mean 
bulk  velocity)  by  comparing  the  results  to  those  of  physical  experiments 
and  previous  numerical  studies.  This  comparison  shows  that  with  a fourth- 
order,  symmetry-preserving  method  a 64  x 64  x 32  grid  suffices  to  perform 
an  accurate  direct  numerical  simulation. 


1.  Introduction 

The  smallest  scales  of  motion  in  a turbulent  flow  result  from  a subtle 
balance  between  convective  transport  and  diffusive  dissipation.  In  math- 
ematical terms,  the  balance  is  an  interplay  between  two  differential  oper- 
ators differing  in  symmetry:  the  convective  derivative  is  skew-symmetric, 
whereas  diffusion  is  governed  by  a symmetric,  definite  operator.  With  this 
in  mind,  we  have  developed  a spatial  discretization  method  which  preserves 
the  symmetries  of  the  balancing  differential  operators.  That  is,  convection 
is  approximated  by  a skew-symmetric  discrete  operator,  and  diffusion  is  dis- 
cretized by  a symmetric,  definite  operator.  Second-order  and  fourth-order 
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versions  have  been  developed  thus  far,  applicable  to  structured  nonuniform 
grids.  The  resulting  semi-discrete  representation  conserves  energy  exactly 
in  the  absence  of  physical  dissipation.  For  finite  Reynolds  numbers,  i.e.  in 
the  presence  of  physical  dissipation,  the  kinetic  energy  of  any  discrete  solu- 
tion decreases  unconditionally  in  time.  Therefore,  a symmetry-preserving, 
spatial  discretization  is  stable  on  any  grid,  and  we  need  not  add  an  artificial 
demping  mechanism  that  will  inevitably  interfere  with  the  subtle  balance 
between  convection  and  diffusion  at  the  smallest  length  scales.  This  forms 
our  motivation  to  investigate  symmetry-preserving  discretizations  for  direct 
numerical  simulation  (DNS)  of  turbulent  flow.  Because  stability  is  not  an 
issue,  the  main  question  becomes  how  accurate  is  a symmetry-preserving 
discretization,  or  stated  otherwise,  how  coarse  may  the  grid  be  for  a DNS? 
We  will  address  this  question  in  Section  2 by  evaluating  the  results  for  a 
turbulent  flow  in  a channel  at  Re=5,600.  We  will  kick  off  by  sketching  the 
main  lines  of  symmetry-preserving  discretization  (Section  1).  For  a more 
detailed  discussion,  we  refer  to  Verstappen  and  Veldman  (1998,  2002).  Con- 
servation properties  of  numerical  schemes  for  the  Navier-Stokes  equations 
are  currently  also  pursued  at  other  research  institutes,  see  e.g.  Hyman  et  al. 
(1992),  Morinishi  et  al.  (1998),  Vasilyev  (2000),  Twerda  (2000)  and  Ducros 
et  al.  (2000). 

2.  Symmetry-preserving  discretization 

The  temporal  evolution  of  the  discrete  velocity  vector  uh  is  governed  by  a 
finite- volume  discretization  of  the  incompressible  Navier-Stokes  equations: 

n^  + C(uh)uh  + Duh-M*ph  = 0 Muh  = 0,  (1) 

where  the  vector  ph  denotes  the  discrete  pressure,  O is  a (positive-definite) 
diagonal  matrix  representing  the  sizes  of  the  control  volumes,  C ( Uh ) is  built 
from  the  convective  flux  contributions  through  the  control  faces,  D contains 
the  diffusive  fluxes,  and  M is  the  coefficient  matrix  of  the  discretization  of 
the  integral  form  of  the  law  of  conservation  of  mass.  The  coefficient  matrices 
C ( Uh ) and  D are  constructed  such  that 

C ( Uh ) + C*  ( Uh ) = 0,  D + D*  positive- definite.  (2) 

The  symmetry  condition  on  the  coefficient  matrix  C ( Uh ) reflects  that 
C (uh)  represents  a discrete  gradient  operator:  its  null  space  consists  of 
the  constant  vectors  (provided  that  the  consistency  condition  C (tt/j)  1=0 
is  satisfied)  and  C ( Uh ) is  skew-symmetric  like  a first-order  differential  op- 
erator. The  coefficient  matrix  D of  the  discrete  diffusive  operator  inherits 
its  positive-definiteness  from  the  underlying  diffusive,  differential  operator 
-V  • V/Re. 
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The  semi-discretization  (l)-(2)  is  conservative  and  stable.  The  total 
mass  is  trivially  conserved,  and  the  same  holds  for  the  total  amount  of 
momentum  (provided  that  the  discretization  is  exact  for  Uh  = 1).  The 
evolution  of  the  discrete  energy  u*hfluh  of  any  solution  Uh  of  (l)-(2)  is 
governed  by 

ft(u*hCluh)  (=]  -uUC  + Cnuh-ul(D  + D*)uh  + 2plM^ 

-0 

= ~u*h(D  + D*)uh  < 0,  (3) 

where  we  have  used  the  skew-symmetry  of  C (it/,).  The  right-hand  side  is 
zero  if  and  only  if  Uh  = 0,  or  D + D*  = 0.  Thus,  the  energy  is  conserved 
if  the  diffusion  is  turned  off.  Note  that  the  pressure  term  M*ph  in  (1) 
does  not  affect  the  evolution  of  the  total  kinetic  energy  (on  condition  that 
Muh  = 0),  because  the  discrete  pressure  gradient  is  represented  by  the 
transpose  of  the  coefficient  matrix  M of  the  incompressibility  constraint. 

With  diffusion  (that  is  for  D/0)  the  right-hand  side  of  (3)  is  negative 
for  all  Uh  ^ 0 provided  that  D + D*  is  positive-definite.  So,  the  energy  of 
the  discrete  system  (1)  decreases  in  time  if  (2)  is  satisfied.  The  semi-discrete 
system  (1)  is  stable  under  this  symmetry  condition:  a solution  of  (1)  can 
be  obtained  on  any  grid,  and  there  is  no  need  to  add  an  artificial  damping 
mechanism  to  stabilize  the  spatial  discretization. 

Since  these  favorable  conservation  and  stability  properties  are  directly 
related  to  the  symmetries  of  the  coefficient  matrices  in  (1),  we  want  to 
construct  these  matrices  such  that  they  fulfil  (2).  To  illustrate  the  way  in 
which  this  may  be  achieved,  we  consider  the  approximation  of  the  first-order 
derivative  in  one  spatial  dimension.  The  traditional  method  maximizes  the 
(formal)  order  of  the  local  truncation  error.  On  a stencil  consisting  of  three 
points,  this  leads  to  the  second-order  approximation 

dxu{xi ) « {r~lUi+ 1 - (r^1  - , (4) 

where 

^i”1  = \ (Si+l  ~ xi- 1)  and  ri  = (*i+l  ~ xi)/(xi  ~ xi- 1)- 

The  essence  of  our  method  is  that  the  first-order  derivative  dxu(xi ) is  ap- 
proximated by  a discrete  operator  f2_1C,  where  the  coefficient  matrix  C 
is  skew-symmetric: 

dxu{Xi)  ~ 2^2  (wi+l  ^i+l)  • (^) 

The  two  ways  of  discretization  are  illustrated  in  Figure  1. 
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Figure  1.  Two  ways  of  approximating  dxu.  In  the  left-hand  figure  the  derivative  is  ap- 
proximated by  means  of  a Lagrangian  interpolation,  that  is  by  Eq.  (4).  In  the  right-hand 
figure  the  symmetry-preserving  discretization  (5)  is  applied. 


As  the  diagonal-entry  of  operator  in  the  right-hand  side  of  (4)  is  non-zero 
for  ri  1,  the  standard  discretization  method  breaks  the  skew- symmetry 
on  a nonuniform  grid.  Consequently,  the  standard  method  does  not  conserve 
the  energy  and  is  not  conditionally  stable  on  nonuniform  meshes. 

The  local  truncation  error  of  the  symmetry-preserving  discretization 

Th{xi)  = 2(<fai+i  - 8xi)dxxu{xi)  + 0(&4,ax),  (6) 

is  first-order,  unless  the  grid  is  (almost)  uniform.  Given  stability,  a sufficient 
condition  for  second-order  accuracy  of  the  discrete  solution  u*  is  that  the 
local  truncation  error  Th  be  second  order.  Then  the  error  e/t  in  the  solution 
Uhi  given  by  0_1Ce/l  = Th,  is  second-order.  Yet,  this  is  not  a necessary  con- 
dition, as  is  emphasized  by  Manteufel  and  White  (1986).  They  have  proven 
that  the  symmetry-preserving  approximation  yields  second-order  accurate 
solutions  on  uniform  as  well  as  on  nonuniform  meshes,  even  though  its  local 
truncation  error  Th  is  formally  only  first-order  on  nonuniform  meshes. 

The  accuracy  of  the  basic  scheme  (5)  may  be  improved  by  means  of 
a Richardson  extrapolation,  just  like  in  Antonopoulos-Domis  (1981).  This 
results  into  the  following,  fourth-order  accurate  discretization: 

dxu(xi)  « (-Ui+2  + 8ui+i  - 8uj_i  + Ui-2) , 


where 


— 2 { 3^+2  4"  8Xi^.\  8Xj_i  T Xj— 2)  . 


The  diffusive  operator  undergoes  a similar  treatment,  leading  to 


Qidxx'U'(Xi)  — 8 


( ^i+1  ui 

X{ 


'U'i  'U'i— 1 

X{—\ 


( 2 'U'i  2 \ 

\*^z-+*2  %i— 2/ 
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Next,  we  will  compare  the  symmetry-preserving  discretization  with  the 
traditional  discretization  methods  based  on  Lagrange  interpolation  (mini- 
mizing local  truncation  error)  for  a steady-state  solution  of  the  convection- 
diffusion  equation 

dtu  + udxu  — dxxu/ Re  = 0. 

Since  on  uniform  grids  the  methods  are  equal,  we  choose  an  example  with  a 
boundary-layer  character,  requiring  grid  refinement  near  the  outflow  bound- 
ary x = 1.  This  is  achieved  by  imposing  the  boundary  conditions  u(0,  t)  = 0 
and  tt(l,  t)  = 1.  The  parameters  are  set  equal  to  u = 1 and  Re  = 1, 000. 

Grid  refinement  has  been  carried  out  on  an  exponentially  stretched  grid, 
with  half  the  grid  points  in  the  thin  boundary  layer  of  thickness  10/Re  near 
x = 1.  Four  discretization  methods  have  been  investigated: 

— The  traditional  Lagrangian  second-order  method  (2L)  and  its  fourth- 
order  counterpart  (denoted  by  4L)  where  we  have  implemented  exact 
boundary  conditions  to  circumvent  the  problem  of  a difference  molecule 
that  is  too  large  near  the  boundary; 

— The  second-order  (2S)  and  fourth-order  (4S)  symmetry-preserving  meth- 
ods. 

We  form  the  vector  wexact  by  restricting  the  analytical  solution  to  the 
grid  points,  and  monitor  the  global  discretization  error  defined  by  | — 

^exactlU  (where  the  norm  is  the  kinetic  energy  norm).  Figure  2 shows  the 
global  error  as  a function  of  the  mean  mesh  size  1/iV,  where  N is  the  number 
of  grid  points. 

A number  of  observations  can  be  made. 

— For  all  grid  sizes  the  Lagrangian  discretization  appears  to  be  less  ac- 
curate than  its  symmetry-preserving  alternative. 

— For  coarser  grids  the  4th-order  Lagrangian  method  is  not  even  as  ac- 
curate as  its  2nd-order  Lagrangian  relative.  Similar  observations  have 
been  made  frequently,  and  this  explains  why  thusfar  4th-order  dis- 
cretization has  not  been  very  popular. 

— The  symmetry-preserving  methods  already  behave  nicely  on  coarse 
grids.  They  display  a regular  monotone  behaviour  upon  grid  refine- 
ment. Moreover,  the  discretization  error  of  4S  (2S)  picks  up  its  final 
slope  at  much  coarser  grids  then  4L  (2L).  As  in  turbulent-flow  simu- 
lations one  will  always  have  to  cope  with  limitations  on  the  affordable 
number  of  grid  points,  methods  that  are  less  sensitive  in  this  respect 
are  preferable. 

— Also  note  that  for  a given  accuracy  (say  10-5),  the  grid  size  of  the 
4th-order  symmetry-preserving  method  can  be  chosen  roughly  three 
times  larger  than  that  of  the  4th-order  Lagrangian  method! 
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number  of  unstable  eigenvalues 


Figure  2.  The  left-hand  figure  shows  the  global  error  as  a function  of  the  mean  mesh  size 
on  an  exponential  grid  with  half  of  the  grid  points  inside  a boundary  layer  of  thickness 
10/Re.  Four  methods  are  shown:  2L  and  4L  (2nd-  and  4th-order  Lagrangian),  2S  and 
4S  (2nd-  and  4th-order  symmetry-preserving).  The  right-hand  figure  depicts  the  number 
of  eigenvalues  of  the  Lagrangian  methods  2L  and  4L  located  in  the  unstable  halfplane. 
Only  the  Lagrangian  methods  are  shown,  since  the  symmetry-preserving  discretization 
keeps  all  the  eigenvalues  in  the  stable  halfplane. 


- The  fourth-order  Lagrangian  method  nearly  breaks  down  for  N = 28 
where  the  stretching  factor  is  0.72  (which  is  not  extreme).  This  is  due 
to  an  eigenvalue  moving  from  the  unstable  half  plane  (for  low  values 
of  N),  towards  the  stable  half  plane  (for  higher  N),  which  crosses  the 
imaginary  axis  close  to  the  origin,  making  the  coefficient  matrix  almost 
singular.  When  one  or  more  eigenvalues  of  the  coefficient  matrix  are 
located  in  the  unstable  halfplane,  the  corresponding  time-dependent, 
semi-discrete  system  is  unstable,  and  can  not  be  integrated  in  the  time 
domain.  In  the  above  examples  we  have  computed  the  discrete  steady- 
state  by  a direct  matrix  solver  to  avoid  this  problem. 

For  details  concerning  the  application  to  the  three-dimensional,  incom- 
presible  Navier-Stokes  equation,  we  refer  to  Verstappen  and  Veldman  (1998, 
2002).  On  a uniform  grid,  the  second  order  scheme  proposed  by  Harlow 
and  Welsh  (1965)  preserves  the  symmetries  of  the  convective  and  diffusive 
operator.  In  outline,  we  have  generalized  Harlow  and  Welsh’s  scheme  to 
nonuniform  meshes  in  such  a manner  that  the  symmetries  are  not  broken, 
and  apply  Richardson  extrapolation  to  improve  the  order  of  accuracy.  A 
variant  of  our  approach  for  collocated  grids  has  been  developed  at  Delft 
University  (Twerda,  2000). 
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3.  A challenging  test-case:  turbulent  channel  flow 


In  this  section,  the  symmetry-preserving  discretization  is  tested  for  turbu- 
lent channel  flow.  The  Reynolds  number  is  set  equal  to  Re  = 5,600  (based  on 
the  channel  width  and  the  bulk  velocity),  a Reynolds  number  at  which  di- 
rect numerical  simulations  have  been  performed  by  several  research  groups; 
see  Kim  et  al.  (1987),  Gilbert  and  Kleiser  (1991),  Kuroda  et  al.  (1995).  In 
addition  we  can  compare  the  numerical  results  to  experimental  data  from 
Kreplin  and  Eckelmann  (1979). 

As  usual,  the  flow  is  assumed  to  be  periodic  in  the  stream-  and  spanwise 
direction.  Consequently,  the  computational  domain  may  be  confined  to  a 
channel  unit  of  dimension  2n  x 1 x n,  where  the  width  of  the  channel  is 
normalized.  All  computations  presented  in  this  section  have  been  performed 
with  64  (uniformly  distributed)  streamwise  grid  points  and  32  (uniformly 
distributed)  spanwise  points.  In  the  lower-half  of  the  channel,  the  wall- 
normal  grid  points  are  computed  according  to 


sinh  jjj/Ny) 
2 sinh  (7/2) 


with  j = 0,  l,...,Ny/2, 


where  Ny  denotes  the  number  of  grid  points  in  the  wall-normal  direction. 
The  stretching  parameter  7 is  taken  equal  to  6.5.  The  grid  points  in  the 
upper-half  are  computed  by  means  of  symmetry. 

The  temporal  integration  of  (1)  is  performed  with  the  help  of  a one- 
leg  method  that  is  tuned  to  improve  its  convective  stability  (Verstap- 
pen  and  Veldman,  1997).  The  non-dimensional  time  step  is  set  equal  to 
St  — 1.25  10-3.  Mean  values  of  computational  results  are  obtained  by  av- 
eraging the  results  over  the  directions  of  periodicity,  the  two  symmetrical 
halves  of  the  channel,  and  over  time.  The  averaging  over  time  starts  after  a 
start-up  period.  The  start-up  period  as  well  as  the  time-span  over  which  the 
results  are  averaged,  1500  non-dimensional  time-units,  are  identical  for  all 
the  results  shown  is  this  section.  Figure  3 shows  a comparison  of  the  mean 
velocity  profile  as  obtained  from  our  fourth-order  symmetry-preserving  sim- 
ulation (Ny  = 64)  with  those  of  other  direct  numerical  simulations.  Here 
it  may  be  stressed  that  the  grids  used  by  the  DNS’s  that  we  compare  with 
have  typically  about  1283  grid  points,  that  is  16  times  more  grid  points 
than  our  grid  has.  Nevertheless,  the  agreement  is  excellent. 

To  investigate  the  convergence  of  the  fourth-order  method  upon  grid 
refinement,  we  have  monitored  the  skin  friction  coefficient  Cf  as  obtained 
from  simulations  on  four  different  grids.  We  will  denote  these  grids  by  A, 
B,  C and  D.  Their  spacings  differ  only  in  the  direction  normal  to  the  wall. 
They  have  Ny  = 96  (grid  A),  Ny  = 64  (B),  Ny  = 56  (C)  and  Ny  = 48  (D) 
points  in  the  wall-normal  direction,  respectively.  The  first  (counted  from 
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Figure  3.  The  mean  streamwise  velocity  u+  versus  y+.  The  dashed  lines  represent  the 
law  of  the  wall  and  the  log  law.  The  markers  represent  DNS-results  that  are  taken  from 
the  ERCOFTAC  Database. 


the  wall)  grid  line  used  for  the  convergence  study  is  located  at  « 0.95 
(grid  A),  yf  » 1.4  (B),  y±  « 1.6  (C),  and  y+  « 1.9  (D),  respectively. 
Figure  4 displays  the  skin  friction  coefficient  Cf  as  function  of  the  fourth 


Figure  4 ■ Convergence  of  the  skin  friction  coefficient  Cf  upon  grid  refinement.  The 
figure  displays  Cf  versus  the  fourth  power  of  the  first  grid  point  yf. 
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power  of  yf.  The  convergence  study  shows  that  the  discretization  scheme  is 
indeed  fourth-order  accurate  (on  a nonuniform  mesh).  This  indicates  that 
the  underlying  physics  is  resolved  when  48  or  more  grid  points  are  used  in 
the  wall  normal  direction.  In  terms  of  the  local  grid  spacing  (measured  by 
yf),  the  skin  friction  coefficient  is  approximately  given  by  Cf  — 0.00836  — 
0.000004(yJf )4.  The  extrapolated  value  at  yf  = 0 lies  in  between  the  Cf 
reported  by  Kim  et  al.  (1987)  (Cf  = 0.00818)  and  Dean’s  correlation  of 
Cf  = 0.073  Re~l /4  = 0.00844  (Dean,  1978). 

The  convergence  of  the  fluctuating  streamwise  velocity  near  the  wall 
(0  < y+  < 20)  is  presented  in  Figure  5.  Here,  we  have  added  results  obtained 
on  three  still  coarser  grids  (with  Ny  = 32,  Ny  = 24  and  Ny  = 16  points 
in  the  wall- normal  direction,  respectively),  since  the  results  on  the  grids 
A,  B,  C and  D fall  almost  on  top  of  each  other.  The  coarsest  grid,  with 
only  Ny  = 16  points  to  cover  the  channel  width,  is  coarser  than  most  of 
the  grids  used  to  perform  a large-eddy  simulation  (LES)  of  this  turbulent 
flow.  Nevertheless,  the  64  x 16  x 32  solution  is  not  that  far  off  the  solution 
on  finer  grids,  in  the  near  wall  region.  Further  away  from  the  wall,  the 
turbulent  fluctuations  predicted  on  the  coarse  grids  (Ny  < 32)  become  too 
high  compared  to  the  fine  grid  solutions,  as  is  shown  in  Figure  6. 

The  solution  on  the  64  x 24  x 32,  for  example,  forms  an  excellent  starting 
point  for  a large-eddy  simulation.  The  root-mean-square  of  the  fluctuating 
streamwise  velocity  is  not  far  of  the  fine  grid  solution,  and  viewed  through 
physical  glasses,  the  energy  of  the  resolved  scales  of  motion,  the  coarse 
grid  (Ny  = 24)  solution,  is  convected  in  a stable  manner,  because  it  is 
conserved  by  the  discrete  convective  operator.  Therefore,  we  think  that  the 
symmetry-preserving  discretization  forms  a solid  basis  for  testing  sub-grid 
scale  models.  The  discrete  convective  operator  transports  energy  from  a 
resolved  scale  of  motion  to  other  resolved  scales  without  dissipating  any 
energy,  as  it  should  do  from  a physical  point  of  view.  The  test  for  a sub- 
grid scale  model  then  reads:  does  the  addition  of  the  dissipative  sub-grid 
model  to  the  conservative  convection  of  the  resolved  scales  reduce  the  error 
in  the  computation  of  itrms. 

The  results  for  the  fluctuating  streamwise  velocity  urms  are  compared 
to  the  experimental  data  of  Kreplin  and  Eckelmann  (1979)  and  to  the 
numerical  data  of  Kim  et  al.  (1987)  in  Fig.  7.  This  comparison  confirms 
that  the  fourth-order,  symmetry-preserving  method  is  more  accurate  than 
the  second-order  method.  With  48  or  more  grid  points  in  the  wall  normal 
direction,  the  root-mean-square  of  the  fluctuating  velocity  obtained  by  the 
fourth-order  method  is  in  close  agreement  with  that  computed  by  Kim  et 
al.  (1987)  for  y+  > 20  (Figure  7 shows  this  only  for  y+  up  to  40;  yet, 
the  agreement  is  also  excellent  for  y+  > 40).  In  the  vicinity  of  the  wall 
(y+  < 20),  the  velocity  fluctuations  of  the  fourth-order  simulation  method 
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Figure  5.  The  root-mean-square  velocity  fluctuations  normalized  by  the  wall  shear 
velocity  as  function  of  the  wall  coordinate  y+  on  various  grids  for  y+  < 20.  The  markers 
correspond  to  the  results  obtained  in  the  grid  points.  The  solution  on  grid  B is  also 
represented  by  a continuous  line. 
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Figure  6.  The  root-mean-square  velocity  fluctuations  normalized  by  the  wall  shear 
velocity  for  y+  < 200  on  various  grids. 


fit  the  experiment  data  nicely,  even  up  to  very  coarse  grids  with  only  24 
grid  points  in  the  wall-normal  direction.  However,  the  turbulence  intensity 


DNS  USING  SYMMETRY-PRESERVING  DISCRETIZATION  217 


Figure  7.  Comparison  of  the  mean-square  of  the  streamwise  fluctuating  velocity  as 
function  of  y+. 


in  the  sublayer  (0  < y+  < 5)  predicted  by  the  simulations  is  higher  than 
that  in  the  experiment.  According  to  the  fourth-order  simulation  the  root- 
mean-square  approaches  the  wall  like  urms  ~ 0.38y+  (Ny  — 64).  The  exact 
value  of  this  slope  is  hard  to  pin-point  experimentally.  Hanratty  et  al.  (1977) 
have  fitted  experimental  data  of  several  investigators,  and  thus  came  to  0.3. 
Most  direct  numerical  simulations  yield  higher  values.  Kim  et  al  (1987) 
and  Gilbert  and  Kleiser  (1991)  have  found  slopes  of  0.3637  and  0.3824 
respectively,  which  is  in  close  agreement  with  the  present  findings. 

So,  in  conclusion,  the  results  of  the  fourth-order  symmetry-preserving 
discretization  agree  better  with  the  available  reference  data  than  those  of  its 
second-order  counterpart,  and  with  the  fourth-order  method  a 64  x 64  x 32 
grid  suffices  to  perform  an  accurate  DNS  of  a turbulent  channel  flow  at 
Re=5,600. 
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